home *** CD-ROM | disk | FTP | other *** search
/ Mac Mania 4 / MacMania 4.toast / / Demo's / Igor Demo Pro / 3 PutContentsIn Igor Pro Folder / Technical Notes / Igor Tech Notes / TN020-A Custom Peak Meas / Template v1.3 Folder / procedure < prev    next >
Text File  |  1993-07-15  |  29KB  |  1,242 lines
  1. | Template Peak Measurement Experiment V1.3
  2. | Changes since V1.2: Added Initialize Template macro
  3. Menu "Macros"
  4.     "Initialize Template"
  5.     m_so
  6.     m_mp
  7.     "-"
  8.     Submenu m_b
  9.         m_bi
  10.         m_ba
  11.         m_bd
  12.         m_bf
  13.         m_s
  14.         m_bs
  15.     End
  16.     Submenu m_f
  17.         m_ii
  18.         m_ia
  19.         m_id
  20.         m_iap
  21.         "-"
  22.         m_ff
  23.         "-"
  24.         m_fr
  25.         m_fs
  26.         m_fh
  27.         m_s
  28.     End
  29.     "-"
  30.     m_cug
  31.     m_cw
  32. End
  33. Macro MakePeaks(py,noise,minx,maxx,log,pts)
  34.     String py=g_w
  35.     Variable noise=mp_noise,log=mp_log,pts=mp_pnts,minx=mp_minx,maxx=mp_maxx
  36.     Prompt py,"Output wave name"
  37.     Prompt noise,"peak wave noise (enoise val)"
  38.     Prompt minx,"starting x value"
  39.     Prompt maxx,"ending x value"
  40.     Prompt log,"x wave increment",popup,"linear;logarithmic"
  41.     Prompt pts,"points in output waves",popup,p_pts
  42.  
  43.     Silent 1;PauseUpdate
  44.     mp_noise=noise;mp_log=log;mp_pnts=pts;mp_minx=minx;mp_maxx=maxx
  45.     String ctrX="W_MakeCentersX",ampY="W_MakeAmpsY",widthX="W_MakeWidthsX",px
  46.     String cw="W_MyPkCoefs",wtmp="W_tmp"
  47.     Variable kk1,kk3,last=strlen(py)-1,n=2^pts,npks
  48.     if(cmpstr("y",py[last])==0)
  49.         last-=1
  50.     endif
  51.     px=py[0,last]+"X"
  52.     Mk(py,n);Mk(px,n);$py=0
  53.     SetScale x,minx,maxx,$py,$px
  54.     if(log==2)
  55.         $px=minx*(maxx/minx)^(p/n)
  56.     else
  57.         $px=x
  58.     endif
  59.     WaveStats/Q $ctrX;npks=V_npnts
  60.     iterate (npks)
  61.         n= MyPeakSizesToCoefs($ctrX[i],$ampY[i],$widthX[i],$cw)    |  output cw, for initial guesses
  62.         $py+=MyPeak($cw,$px[p])
  63.     loop
  64.     $py+=enoise(noise)
  65.     g_w=py
  66.     if(log==2)
  67.         g_wx=px
  68.         AppWv(py,px,"")
  69.     else
  70.         g_wx=calc
  71.         AppWv(py,"","")
  72.     endif
  73. EndMacro
  74.  
  75. Macro InitializeTemplate(wty,wtx,norm,type)
  76.     String wty=it_y,wtx=it_x,type=S_PkT
  77.     Variable norm=it_n
  78.     Prompt wty,"template y wave",popup,"cos**2;"+WaveList("*",";","")
  79.     Prompt wtx,"template x wave",popup, none+WaveList("*",";","")
  80.     Prompt norm,"normalize template to amplitude = 1.0?",popup,"yes;no"
  81.     Prompt type,"peak description"
  82.     
  83.     Silent 1;PauseUpdate
  84.     it_y=wty;it_x=wtx;it_n=norm;S_PkT=type    
  85.  
  86.     | Build XY Template
  87.     if(cmpstr(wty,"cos**2")==0)
  88.         it_x="_none_"
  89.         Mk("W_PkT",101);Mk("W_PkTx",101)
  90.         SetScale/I x,-pi,pi,W_PkT,W_PkTx
  91.         W_PkTx=x;W_PkT=1+cos(W_PkTx[p]);
  92.     else
  93.         if(exists(wtx)==1)
  94.             ChkLen(wty,wtx)
  95.             Dup(wty,"W_PkT");Dup(wtx,"W_PkTx")
  96.         else
  97.             Dup(wty,"W_PkT")
  98.             Dup(wty,"W_PkTx");W_PkTx=x
  99.         endif
  100.     endif
  101.     WaveStats/Q W_PkT
  102.     if(norm==1)
  103.         V_Tamp=1;W_PkT/=V_max
  104.     else
  105.         V_Tamp=V_max
  106.     endif
  107.     V_TmaxX=W_PkTx[x2pnt(W_PkT,V_maxloc)]
  108.     Printf "Setting x=0 at template maxima (where x was %g)\r",V_TmaxX
  109.     W_PkTx -= V_TmaxX;V_TmaxX=0
  110.     V_Tarea=AreaXYTPt(W_PkTx,W_PkT,0,numpnts(W_PkT)-1)
  111.     V_Tfwhm=FwhmTmpl()
  112. End
  113.  
  114. Macro InitializeMostEverything(w,wx,wb)
  115.     String w=g_w,wx=g_wx,wb=g_b
  116.     Prompt w,p_w,popup, WaveList("*",";","")
  117.     Prompt wx,p_wx,popup,calc+WaveList("*X", ";","")
  118.     Prompt wb,p_b,popup, none+WaveList("*base*",";","")
  119.  
  120.     Silent 1;PauseUpdate
  121.     g_w=w;g_wx=wx; SBs(wb)
  122.     Redimension/D $w
  123.     SameLen(w,breg);$breg=NaN
  124.     Mk(ampsY,2);Mk(ctrsX,2);Mk(widsX,2);Mk(ctrsP,2);Mk(edgsP,4);
  125.     Mk(areaNB,2);Mk(areaX1,2);Mk(areaX2,2);
  126.     Mk(bdcw,2);Mk(pdcw,2)
  127.     Mk("W_MyPkCoefs",3) | MyPeak cw[0,2]
  128.     InitializeTemplate(it_y,it_x,it_n,S_PkT)
  129.     if(exists(wx)==1)
  130.         Redimension/D$wx
  131.     endif
  132.     if(exists(wb)==1)
  133.         Redimension/D $wb
  134.     endif
  135.     DoWindow/F PeakFitGraph
  136.     AppWv(w,wx,"");HideFittedPeaks()
  137. End
  138.  
  139. Macro InitBaseLineFit(w,wx)
  140.     String w=g_w,wx=g_wx
  141.     Prompt w,p_w,popup, WaveList("*",";","")
  142.     Prompt wx,p_wx,popup,calc+WaveList("*X", ";","")
  143.     
  144.     Silent 1;PauseUpdate
  145.     g_w=w;g_wx=wx
  146.     SameLen(w,breg);AppWv(breg,wx,"");Modify mode($breg)=1;$breg=NaN;AppCrs(w)
  147. End
  148. Macro AddRegionToFit()
  149.     Silent 1;PauseUpdate
  150.     CheckTwoCursors(g_w);$breg[V_start,V_theEnd]=$g_w[p]
  151. End
  152. Macro DeleteRegionFromFit()
  153.     Silent 1;PauseUpdate
  154.     CheckTwoCursors(g_w);$breg[V_start,V_theEnd]=NaN
  155. End
  156.  
  157. Macro FitBaselineAtRegions(w,wx,wr,fit)
  158.     String w=g_w,wx=g_wx,wr=fbar_wr,fit=fbar_fit
  159.     Prompt w,p_w,popup, WaveList("*",";","")
  160.     Prompt wx,p_wx,popup,calc+WaveList("*X", ";","")
  161.     Prompt wr,"region wave (baseline weighting wave)",popup, none+breg
  162.     Prompt fit,"baseline function",popup,S_funcs
  163.     
  164.     Silent 1;PauseUpdate
  165.     g_w=w;g_wx=wx;fbar_wr=wr;fbar_fit=fit
  166.     ChkLen(w,wx);ChkLen(w,wr)
  167.     String ow="W_BaselineFit",wtmp="W_tmp",opts="";SameLen(w,ow)
  168.     Variable t=0,typ=floor(strsearch(S_funcs,fit,0)/10)
  169.     
  170.     if(exists(wx)==1)
  171.         opts+="/X="+wx
  172.     endif
  173.     if(exists(wr)==1)
  174.         Dup(wr,wtmp)
  175.         $wtmp=(numtype($wtmp)==0)
  176.         opts+="/W="+wtmp
  177.         ar_ex=3
  178.     else
  179.         ar_ex=1
  180.     endif
  181.     t=str2num(fit[7,8])-1
  182.     String command="CurveFit "+fit[1,7]+", $w"+opts
  183.     Execute command
  184.     KillWv(wtmp)
  185.     Mk(bcw,2);$bcw={NaN,K0,K1,K2,K3,K4};Redimension/N=(2+t) $bcw
  186.     Mk(bdcw,2);$bdcw={1,0,numpnts($w)-1,typ,t};Dup(bdcw,"W_PM")
  187.     if(exists(wx)==1)
  188.         $ow=PolyMorph($bcw,$wx[p])
  189.     else
  190.         $ow=PolyMorph($bcw,x)
  191.     endif
  192.     AppWv(ow,wx,"")
  193.     SBs(ow)
  194.     ar_wfit=ow
  195. End
  196.  
  197. Macro SubtractBaseline(w,wx,wb,ow)
  198.     String w=g_w,wx=g_wx,wb=rb_b,ow=rb_ow
  199.     Prompt w,p_w,popup, WaveList("*",";","")
  200.     Prompt wx,p_wx,popup,calc+WaveList("*X", ";","")
  201.     Prompt wb,p_b,popup, "_From Fit Peaks_;"+WaveList("*base*",";","")
  202.     Prompt ow,"output baseline-corrected wave",popup,new+WaveList("*",";","")
  203.  
  204.     Silent 1;PauseUpdate
  205.     g_w=w;g_wx=wx
  206.     ChkLen(w,wx);ChkLen(w,wb)
  207.     String cw,dcw,tcw="W_tmp0",tdcw="W_PM"
  208.     if(exists(ow)!=1)
  209.         NewWv(w,"NoBase");ow=S_Wave
  210.     else
  211.         SameLen(w,ow)
  212.     endif
  213.     rb_ow=ow
  214.     if(exists(wb)!=1)
  215.         cw=bcw
  216.         dcw=bdcw
  217.         if(numtype($dcw[1])!=0)
  218.             Abort "Run Fit Peaks first!"
  219.         endif
  220.         Variable terms=$dcw[4],s=$dcw[1],e=$dcw[2],funcs=$dcw[0]
  221.         Dup(cw,tcw);Redimension/N=(2+terms) $cw
  222.         Dup(dcw,tdcw);$tdcw[0]=1
  223.         wb="W_tmp";SameLen(w,wb);$wb=NaN
  224.         if(exists(wx)==1)
  225.             $wb[s,e]=PolyMorph($tcw,$wx[p])
  226.         else
  227.             $wb[s,e]=PolyMorph($tcw,x)
  228.         endif
  229.     endif
  230.  
  231.     $ow=$w-$wb
  232.     AppWv(ow,wx,"");Modify zero(left)=1
  233.     KillWv(wb);KillWv(tcw)
  234.     g_w=ow
  235.     SBs("_None_")
  236. End
  237.  
  238. Macro InitIdentifyPeaks(w,wx,wb,pol,box)
  239.     String w=g_w,wx=g_wx,wb=g_b
  240.     Variable box=g_bx,pol=mpr_pol
  241.     Prompt w,p_w,popup, WaveList("*",";","")
  242.     Prompt wx,p_wx,popup,calc+WaveList("*X", ";","")
  243.     Prompt wb,p_b,popup, none+WaveList("*base*",";","")
  244.     Prompt pol,"peak polarity", popup,"positive;negative"
  245.     Prompt box,p_bx
  246.     
  247.     Silent 1;PauseUpdate
  248.     g_w=w;g_wx=wx;SBs(wb)
  249.     box=abs(box);mpr_pol=pol;g_bx=box
  250.  
  251.     ChkLen(w,wx);ChkLen(w,wb)
  252.     Mk(ampsY,2);Mk(ctrsX,2);Mk(ctrsP,2);Mk(widsX,2);Mk(edgsP,4)
  253.     DoWindow/F PeakFitGraph
  254.     if(V_Flag==0)
  255.         AppWv(w,wx,"");DoWindow/C PeakFitGraph
  256.     endif
  257.     AppWv(w,wx,"");AppWv(wb,wx,"")
  258.     ShowEdges();AppWv(ampsY,ctrsX,"");Modify mode($ampsY)=8,marker($ampsY)=18
  259.     AppCrs(w);HideFittedPeaks()
  260. End
  261.  
  262. Macro IdentifyPeakWithCursors()
  263.     Silent 1;PauseUpdate
  264.     String w=g_w,wx=g_wx,wb=g_b
  265.     CheckTwoCursors(w)
  266.     Variable s=V_start,en=V_theEnd,npks,ctr,box=g_bx,pol=mpr_pol
  267.     String wtmp="W_tmp",wtmp2="W_tmp2",e0="W_tmp0",e1="W_tmp1"
  268.     SplitEdges(e0,e1);npks=V_npnts
  269.     if(exists(wx)==1)
  270.         box *= sign($wx[en]-$wx[s])
  271.         FindPX(w,wx,box+$wx[s]);box=V_peakP
  272.         FindPX(w,wx,$wx[s]);box=abs(V_peakP-box)
  273.     else
  274.         box*=sign(rightx($w)-leftx($w))
  275.         box=x2pnt($w,box+leftx($w))
  276.     endif
  277.     box=max(1,box);Print "Averaging ",box," points..."
  278.     Duplicate/O/R=[s,en] $w,$wtmp;$wtmp=MyBoxSmooth($w,p+s,box)
  279.     WaveStats/Q $wtmp | Find peak's center pt
  280.     if(pol==1)
  281.         ctr=V_maxLoc
  282.     else
  283.         ctr=V_minLoc
  284.     endif
  285.     ctr=x2pnt($wtmp,ctr)+s
  286.     InsertPoints 0,1,$ctrsP,$ctrsX,$ampsY,$widsX,$e0,$e1
  287.     InsertPoints 0,2,$edgsP
  288.     $ctrsP[0]=ctr
  289.     $e0[0]=s;$e1[0]=en
  290.     $ampsY[0]=$wtmp[ctr-s]
  291.     KillWv(wtmp);KillWv(wtmp2)
  292.     if(exists(wb)==1)
  293.         $ampsY[0]-=MyBoxSmooth($wb,ctr,box)
  294.     endif
  295.     if(exists(wx)==1)
  296.         $ctrsX[0]=$wx[ctr]
  297.     else
  298.         $ctrsX[0]=pnt2x($w,ctr)
  299.     endif
  300.     $widsX[0]=abs(V_theEndX-V_startX)
  301.     MergeEdges(e0,e1);ShowEdges()
  302. End
  303.  
  304. Macro DeletePksBetweenCursors()
  305.     Silent 1;PauseUpdate
  306.     String w=g_w,wx=g_wx,wb=g_b
  307.     CheckTwoCursors(w)
  308.     Variable s=V_start,en=V_theEnd,npks,ctr,box=g_bx,pol=mpr_pol
  309.     String wtmp="W_tmp",wtmp2="W_tmp2",e0="W_tmp0",e1="W_tmp1"
  310.     SplitEdges(e0,e1);npks=V_npnts
  311.     iterate (npks)
  312.         if(($ctrsP[i] >= s)%&($ctrsP[i] <= en))
  313.             $ctrsP[i]=NaN;$ctrsX[i]=NaN;$ampsY[i]=NaN;$widsX[i]=NaN
  314.             $e0[i]=NaN;$e1[i]=NaN| MyPeak        
  315.         endif
  316.     loop
  317.     MergeEdges(e0,e1);ShowEdges()
  318.     if(V_npnts<1)
  319.         DoAlert 0, "no peaks left"
  320.     endif
  321. End
  322.  
  323. | Find peaks by analyzing peak data's smoothed derivatives (WARNING: sensitive to noise)
  324. Macro AutoIdentifyPeaks(what,w,wx,wb,mina,pol,box,extent,maxWidth)
  325.     String w=g_w,wx=g_wx,wb=g_b
  326.     Variable mina=fp_minamp
  327.     Variable what=tfp_what,pol=tfp_pol,box=g_bx,extent=tfp_extent,maxWidth=tfp_mw
  328.     Prompt what,"do what with found peaks?",popup, p_aipw
  329.     Prompt w,p_w,popup,WaveList("*", ";","")
  330.     Prompt wx,p_wx,popup,calc+WaveList("*X", ";","")
  331.     Prompt wb,p_b,popup,none+WaveList("*base*", ";","")
  332.     Prompt mina,"min pk amplitude, or 0 for auto-level"
  333.     Prompt pol,"peak polarity", popup,"positive;negative" 
  334.     Prompt box,p_bx
  335.     Prompt extent,"extent of Y wave to search",popup,"entire wave;between cursors"
  336.     Prompt maxWidth,"maximum peak x width (INF for no max)"
  337.  
  338.     Silent 1;PauseUpdate
  339.     tfp_what=what;g_w=w;g_wx=wx;SBs(wb)
  340.     box=abs(box)
  341.     tfp_pol=pol;g_bx=box;tfp_extent=extent
  342.     String cmd="FindAPeak",lvlw="W_ShowMinPkLvl"
  343.     ChkLen(w,wx);ChkLen(w,wb)
  344.     Variable s=0,en=numpnts($w)-1
  345.     if(extent==2) | use cursors
  346.         CheckTwoCursors(w)
  347.         s=V_start
  348.         en=V_theEnd
  349.     endif
  350.     if(exists(wx)==1)
  351.         box *= sign($wx[en]-$wx[s])
  352.         FindPX(w,wx,box+$wx[s]);box=V_peakP
  353.         FindPX(w,wx,$wx[s]);box=abs(V_peakP-box)
  354.     else
  355.         box*=sign(rightx($w)-leftx($w))
  356.         box=x2pnt($w,box+leftx($w))
  357.     endif
  358.     box=max(1,box);Print "Averaging ",box," points..."
  359.     if(mina==0)
  360.         Dup(w,lvlw)
  361.         if(exists(wb)==1)
  362.             $lvlw-=$wb
  363.         endif
  364.         Smooth 3, $lvlw
  365.         WaveStats/Q/R=[s,en] $lvlw
  366.         mina=V_sdev*0.5 |Heuristic
  367.         if (pol==2)
  368.             mina *= -1
  369.         endif
  370.         Print "Auto min peak amplitude = ",mina
  371.     endif
  372.     fp_minamp=mina
  373.     SameLen(w,lvlw);$lvlw=NaN
  374.     if(exists(wb)==1)
  375.         cmd+="/B="+wb
  376.         $lvlw[s,en]=mina+$wb[p]
  377.     else
  378.         $lvlw[s,en]=mina
  379.     endif
  380.     AppWv(lvlw,wx,"");ResumeUpdate
  381.     PauseUpdate
  382.     cmd+=" mina,pol,box,$w[rt,lf]"    |r=>l search
  383.     String e0="W_tmp0",e1="W_tmp1"
  384.  
  385.     if(what==2)
  386.         Mk(ampsY,2);Mk(ctrsX,2);Mk(ctrsP,2);Mk(widsX,2);Mk(edgsP,4) | MyPeak
  387.     endif
  388.  
  389.     Variable lf=s,rt=en,fl,ce,dX,npks=0
  390.     do
  391.         Execute cmd        | FindAPeak...
  392.         if(V_Flag==0)
  393.             | Insert a new Peak
  394.             InsertPoints 0,1,$ctrsP,$ctrsX,$ampsY,$widsX |MyPeak
  395.             InsertPoints 0,2,$edgsP
  396.             $ctrsP[0]=V_peakP
  397.             fl=floor(V_peakP);ce=ceil(V_peakP)
  398.             if(exists(wx)==1)    | convert f.p. V_peakP f.p. to wx coordinate
  399.                 if(fl != ce)
  400.                     dX=$wx[ce]-$wx[fl]
  401.                     $ctrsX[0]=$wx[fl]+(V_peakP-fl)*dX
  402.                 else
  403.                     $ctrsX[0]=$wx[V_peakP]
  404.                 endif
  405.             else
  406.                 $ctrsX[0]=V_peakX
  407.             endif
  408.             rt=fl-floor((box+1)/2)
  409.             npks+=1
  410.             Printf "%d...",npks
  411.         endif
  412.     while ((V_Flag==0) %& ((rt-lf)>1))
  413.     AppWv(ampsY,ctrsX,"");Modify mode($ampsY)=8,marker($ampsY)=18;ResumeUpdate
  414.     PauseUpdate
  415.     Printf "\r%d Peaks found...",npks
  416.     if(npks>0)
  417.         Printf "Estimating their sizes..."
  418.         cmd="EstimatePeakSizes"
  419.         if(exists(wx)==1)
  420.             cmd+="/X="+wx
  421.         endif
  422.         if(exists(wb)==1)
  423.             cmd+="/B="+wb
  424.         endif
  425.         cmd+="/E=$edgsP 50,maxWidth,box,npks,$ctrsP,$w,$ampsY,$widsX"
  426.         Execute cmd
  427.     endif
  428.     SplitEdges(e0,e1);MergeEdges(e0,e1)
  429.     Print "Done"
  430.     ShowEdges()
  431. End
  432.  
  433. Proc SplitEdges(e0,e1)
  434.     String e0,e1
  435.     Variable pts=numpnts($ctrsX)
  436.     Make/O/D/N=(pts) $e0,$e1;$e0=$edgsP[0+p*2];$e1=$edgsP[1+p*2]
  437.     WaveStats/Q $ctrsX | Count non-Nan peaks
  438. End
  439.  
  440. Proc MergeEdges(e0,e1)
  441.     String e0,e1
  442.     Sort $ctrsP,$ctrsX,$ampsY,$widsX,$e0,$e1;Sort $ctrsP,$ctrsP |MyPeak    | NaN's sorted to end of array
  443.     $edgsP[0,;2]=$e0[p/2]
  444.     $edgsP[1,;2]=$e1[(p-1)/2]
  445.     KillWaves $e0,$e1
  446.     WaveStats/Q $ctrsX
  447.     Redimension/N=(V_npnts+2) $ctrsP,$ctrsX,$ampsY,$widsX | MyPeak
  448.     Redimension/N=((V_npnts+2)*2) $edgsP
  449. End
  450.  
  451. Proc ShowEdges()
  452.     Variable b,x0,x1,p0,p1
  453.     String sx="W_ShowEstWidthsX",sy="W_ShowEstWidthsY"
  454.     WaveStats/Q $ampsY
  455.     Mk(sx,V_npnts*5+2);Mk(sy,V_npnts*5+2)
  456.     iterate (V_npnts)
  457.         b=i*5;p0=$edgsP[2*i];p1=$edgsP[2*i+1]
  458.         if(exists(g_wx)==1)
  459.             x0 =$g_wx[p0]
  460.             x1=$g_wx[p1]
  461.         else
  462.             x0=pnt2x($g_w,p0)
  463.             x1=pnt2x($g_w,p1)
  464.         endif
  465.         $sx[b,b+1]=x0
  466.         $sx[b+2,b+3]=x1
  467.         $sy[b,b+3]=0
  468.         $sy[b+1,b+2]=$ampsY[i]
  469.         if(exists(g_b)==1)
  470.             $sy[b,b+1]+=$g_b[p0]
  471.             $sy[b+2,b+3]+=$g_b[p1]
  472.         endif
  473.     loop
  474.     AppWv(sy,sx,"")
  475. End
  476.  
  477. Function/D MyBoxSmooth(w,pp,box)
  478.     Wave/D w
  479.     Variable pp,box | box should be odd
  480.     
  481.     Variable/D result=0
  482.     Variable  hp,top=numpnts(w)-1,terms
  483.     box = 2*trunc(box/2)+1
  484.     pp-=trunc(box/2)
  485.     hp=limit(0,pp+box-1,top)
  486.     pp=limit(0,pp,top)
  487.     terms=hp-pp+1
  488.     do
  489.         result+=w[pp]
  490.         pp+=1
  491.     while (pp<=hp)
  492.     return result/terms
  493. End
  494.  
  495. | W_PM[0]=n, number of functions to fit
  496. | W_PM[1]=start p of peak fit
  497. | W_PM[2]=end p of peak fit
  498. | W_PM[3]=tp, type of baseline
  499. | W_PM[4]=nw, baseline terms, could be 0
  500. | W_PM[5]=tp, type of 1st peak function
  501. |...
  502. | w[0]=NaN
  503. | w[1]=K0
  504. | w[2]=1st  coefficient...
  505. Function/D PolyMorph(w,xx)
  506.     Wave/D w
  507.     Variable/D xx
  508.     
  509.     Variable/D tp,nw, bs=3,n=W_PM[0],st=2,ii,s,r=w[1] | K0
  510.      if(n==0)
  511.         return r
  512.     endif
  513.     do    
  514.         tp=W_PM[bs]
  515.         nw=W_PM[bs+1]
  516.         if(nw>0)
  517.             if(tp<4) | line, poly 3 - poly 5 K1*x+K2*x^2...
  518.                 ii=nw-1
  519.                 s=0
  520.                 do
  521.                     s=w[ii+st]+xx*s
  522.                     ii-=1
  523.                 while (ii >= 0)
  524.                 r+= s*xx
  525.             else
  526.                 if(tp==7) | MyPeak
  527.                     ii=0
  528.                     do
  529.                         W_MyPkCoefs[ii]=w[ii+st]
  530.                         ii+=1
  531.                     while (ii < nw)
  532.                     r+=MyPeak(W_MyPkCoefs,xx)
  533.                 else
  534.                 if(tp==6) | exp
  535.                     r+= w[st]*exp(-w[st+1]*xx)
  536.                 else
  537.                 if(tp==5)    | dblexp
  538.                     r+= w[st]*exp(-w[st+1]*xx)*exp(-w[st+2]*xx)
  539.                 else    | 4, sin
  540.                     r+= w[st]*sin(w[st+1]*xx+w[st+2])
  541.                 endif
  542.                 endif
  543.                 endif
  544.             endif
  545.             st+=nw    
  546.         endif
  547.         bs+=2
  548.         n-=1
  549.     while (n>0)
  550.     return r
  551. End
  552.  
  553. Function MyCoefsToPeakSizes(cw)    | input cw, global outputs |MyPeak
  554.     Wave/D cw
  555.     V_ctr=cw[2]
  556.     V_amp=cw[0]*V_Tamp
  557.     V_fwhm=V_Tfwhm/cw[1] | exact
  558.     V_area=cw[0]/cw[1]*V_TArea
  559.     return 0
  560. End
  561.  
  562. Function MyPeakSizesToCoefs(ctr,amp,fwhm,cw) |output cw, for initial guesses
  563.     Variable/D ctr,amp,fwhm
  564.     Wave/D cw
  565.     cw[0]=amp/V_Tamp
  566.     cw[1]=V_Tfwhm/fwhm | exact
  567.     cw[2]=ctr
  568.     return 3    | #of cw terms
  569. End
  570.  
  571. Function/D MyPeak(w,x) | Template, 3 terms
  572.     Wave/D w
  573.     Variable/D x
  574.     x=limit(w[1]*(x-w[2]),W_PkTx[0],W_PkTx[numpnts(W_PkTx)-1])
  575.     return w[0]*interp(x,W_PkTx,W_PkT)
  576. End
  577.  
  578. | assumes W_PkT has 0 baseline
  579. Function FwhmTmpl()
  580.     Variable n,px0,ylvl,lx,rx
  581.     n=numpnts(W_PkT)
  582.     px0=BinaryLvlSearch(W_PkTx,0,n-1,V_TmaxX)
  583.     ylvl=V_Tamp/2
  584.     lx= XEdgeXY(W_PkTx,W_PkT,0,px0,ylvl)
  585.     rx= XEdgeXY(W_PkTx,W_PkT,px0,n-1,ylvl)    
  586.     return rx-lx
  587. End
  588. | returns x coordinate where y passes through ywave
  589. | assumes: xwave monotonic increasing
  590. Function XEdgeXY(wx,wy,pst,pen,y)
  591.     Wave wx,wy
  592.     Variable pst,pen,y
  593.     Variable pt,lx,y1,y2
  594.     pt= BinaryLvlSearch(wy,pst,pen,y)
  595.     if(pt==-1)
  596.         lx= wx[pst]
  597.     else
  598.     if(pt==-2)
  599.         lx=wx[pen]
  600.     else
  601.         y1=wy[pt]
  602.         y2=wy[min(pt+1,numpnts(wy)-1)]
  603.         if(y1==y2)
  604.             lx= wx[pt]
  605.         else
  606.             lx=wx[pt]+(y-y1)*(wx[pt+1]-wx[pt])/(y2-y1)
  607.         endif
  608.     endif
  609.     endif
  610.     return lx
  611. End
  612.  
  613. | returns point coordinate before or at point where y passes through ywave,
  614. |  -1 if offscale at zero end or -2 if offscale at other end
  615. |
  616. Function BinaryLvlSearch(ywave,pst,pen,y)
  617.     Wave ywave
  618.     Variable pst,pen,y
  619. ;
  620.     variable n= min(1+pen-pst,numpnts(ywave)),tmp
  621.     variable increasing= ywave[pen] > ywave[pst]
  622.  
  623.     variable jl= pst-1, ju= pen+1    | lower and upper bounds
  624.     variable jm                    | midpoint
  625.  
  626.     if( y == ywave[pst] )
  627.         return pst                | special case for endpoints
  628.     endif
  629.     if( y == ywave[pen] )
  630.         return pen
  631.     endif
  632.     
  633.     if( increasing )
  634.         if( y > ywave[pen] )
  635.             return -2
  636.         endif
  637.     else
  638.         if( y < ywave[pen] )
  639.             return -2
  640.         endif
  641.     endif
  642.     do
  643.         if( ju-jl <= 1 )
  644.             break;
  645.         endif
  646.         jm= floor((ju+jl)/2)
  647.         if( (y >= ywave[jm]) == increasing )
  648.             jl= jm
  649.         else
  650.             ju= jm
  651.         endif
  652.     while(1)
  653.     return jl
  654. end
  655. Function AreaXYTPt(wx,wy,pa,pb) | trapezoidal
  656.     Wave wx,wy    | monotonic wx!
  657.     Variable pa,pb    |pa<pb
  658.  
  659.     Variable a=0
  660.     do
  661.         a+=(wy[pa]+wy[pa+1])*(wx[pa+1]-wx[pa])/2
  662.         pa+=1
  663.     while(pa<pb)
  664.     return a
  665. End
  666.  
  667. | This fits peak equations and a baseline equation or wave to the data.
  668. | If the baseline is removed, a constant baseline + peaks will be fit to the data.
  669. Macro FitPeaksAndBaseline(w,wx,extent,wb,wts)
  670.     String w=g_w,wx=g_wx,wb=fpks_b,wts=fpks_weights
  671.     Variable extent=fpks_extent
  672.     Prompt w,p_w,popup, WaveList("*",";","")
  673.     Prompt wx,p_wx,popup,calc+WaveList("*X", ";","")
  674.     Prompt extent,"extent of peak wave to fit",popup,"entire wave, all peaks;wave & peaks within cursors"
  675.     Prompt wb,p_b,popup, none+S_funcs+"_From Fit Baseline_;"+WaveList("*base*",";","")
  676.     Prompt wts,"weighting wave",popup, none+WaveList("*",";","")
  677.     
  678.     Silent 1;PauseUpdate
  679.     g_w=w;g_wx=wx;SBs(wb)
  680.     fpks_extent=extent;fpks_weights=wts
  681.     
  682.     ChkLen(w,wx);ChkLen(w,wb);ChkLen(w,wts)
  683.     ChkLen(ctrsX,ampsY);ChkLen(ctrsX,widsX)
  684.     Variable s=0,en=numpnts($w)-1
  685.     if(extent==2)
  686.         CheckTwoCursors(w)
  687.         s=V_start
  688.         en=V_theEnd
  689.     endif
  690.  
  691.     String tmpw="W_tmp",tw0="W_tmp0",ow="W_PeakFit"
  692.     Mk(pcw,5);Mk(pdcw,5)
  693.     Dup(w,ow);$ow=NaN
  694.     Variable terms=0,pkTyp=0,lim=1,baseTyp = floor(strsearch(S_funcs,wb,0)/10)
  695.     String cmd,opts=" /D="+ow
  696.     WaveStats/Q/R=[s,en] $w
  697.     Variable/D bsgs=V_avg,xdiff=pnt2x($w,en)-pnt2x($w,s)
  698.     if(bsgs==0)
  699.         bsgs=V_tol
  700.     endif
  701.     if(exists(wx)==1)
  702.         opts+=" /X="+wx
  703.         xdiff=$wx[en]-$wx[s]
  704.     endif
  705.     if(exists(wts)==1)
  706.         opts+="/W="+wts
  707.     endif
  708.     AppWv(ow,wx,"")
  709.     | initial baseline guesses
  710.     String hold
  711.     if(baseTyp>=0)| function selected (expect "_poly   1",  "_dblexp 5", etc)
  712.         terms=str2num(wb[7,8])-1 | exclude K0
  713.         cmd="CurveFit/Q/O "+wb[1,7]+", "+w+"[s,en] "+opts
  714.         Execute cmd
  715.         $pcw={NaN,K0,K1,K2,K3,K4}
  716.         $pdcw={1,s,en,baseTyp,terms}
  717.         hold="100000"[0,terms+1]
  718.         afp_b="_From Fit Peaks_"
  719.     else
  720.         if(exists(wb)==1)    | remove base before fit
  721.             Dup(w,tmpw);w=tmpw;$w-=$wb
  722.             terms=0;$pcw={NaN,0};$pdcw={1,s,en,0,0}
  723.             hold="11";bsgs=0 | K0 not varied
  724.         else
  725.             afp_b="_From Fit Peaks_"
  726.             if(cmpstr(wb,"_From Fit Baseline_")==0)    | Use base fit to start
  727.                 Dup(bcw,pcw);Dup(bdcw,pdcw)
  728.                 terms=$pdcw[4]
  729.                 $pdcw[1]=s;$pdcw[2]=en
  730.                 hold="100000"[0,terms+1]
  731.             else | _None_
  732.                 terms=0
  733.                 $pcw={NaN,bsgs}
  734.                 $pdcw={1,s,en,0,0}
  735.                 hold="10" |  K0 varied
  736.             endif
  737.         endif
  738.     endif
  739.     | prevent singular matrix w/ poly base
  740.     if($pdcw[3]<4)
  741.         $pcw[1,1+terms]=$pcw[p]*($pcw[p]!=0)+($pcw[p]==0)*bsgs/(xdiff^(p-1))
  742.     endif
  743.     WaveStats/Q $ctrsX
  744.     Variable npks=V_npnts
  745.     if(npks<1)
  746.         Abort "no peaks in peak center wave "+ctrsX
  747.     endif
  748.     Variable/D st=2+terms,dst,sz,dsz=2,pks=0,n
  749.     String/G holdTerms
  750.     Redimension/N=(5+npks*2) $pdcw
  751.     pkTyp= 7
  752.     sz=numpnts(W_MyPkCoefs)|MyPeak
  753.     holdTerms="00000000"[0,sz-1]
  754.     iterate (npks)
  755.         if(($ctrsP[i]>=s) %& ($ctrsP[i]<=en))
  756.             dst=3+dsz*$pdcw[0]
  757.             Redimension/N=(st+sz) $pcw
  758.             Mk(tw0,sz)
  759.             n= MyPeakSizesToCoefs($ctrsX[i],$ampsY[i],$widsX[i],$tw0)    |  output $tw0, for initial guesses
  760.             $pcw[st,st+sz-1]=$tw0[p-st]
  761.             $pdcw[dst]=pkTyp
  762.             $pdcw[dst+1]=sz
  763.             hold+= holdTerms
  764.             st+=sz
  765.             $pdcw[0]+=1
  766.             pks+=1
  767.         endif
  768.     loop
  769.     if(pks<1)
  770.         Abort "no peaks between cursors"
  771.     else
  772.         Print "Fitting ",pks," peaks"
  773.     endif
  774.     | Fit peaks using initial guesses for coefficients
  775.     Dup(pdcw,"W_PM")
  776.     cmd="FuncFit/H=hold  PolyMorph $pcw,$w[s,en]"+opts
  777.     ResumeUpdate
  778.     Execute cmd        | FuncFit… ("singular matrix…" usually means normalize x values!)
  779.     KillWv(tmpw)
  780.     if(exists(wb)==1)
  781.         $ow[s,en]+=$wb[p]
  782.     endif
  783.     ar_wfit=ow
  784.     ar_ex=extent
  785.     g_keep=w
  786. End
  787.  
  788. Macro Report(title,srt,order,basInf)
  789.     String title=pfr_title, srt=pfr_sort,
  790.     Variable order=pfr_order,basInf=pfr_bi
  791.     Prompt title,"report title"
  792.     Prompt srt,"sort report by",popup,"center;amplitude;width;area"
  793.     Prompt order,"sort order",popup,"ascending;descending"
  794.     Prompt basInf,"baseline info",popup,"include;omit"
  795.     
  796.     Silent 1
  797.     pfr_title=title;pfr_sort=srt;pfr_order=order;pfr_bi=basInf
  798.     String center="W_PkCenters",amplitude="W_PkAmps",width="W_PkFWHMs",area="W_PkAreas",tw2="W_tmp2"
  799.  
  800.     Variable terms=$pdcw[4],btyp=$pdcw[3] | baseline terms, not counting K0
  801.     Variable sz=3,dsz=2,st=2+terms,dst=5,n
  802.     Variable npks=$pdcw[0]-1,np1=npks-1,nrows=max(2,npks)
  803.     String text=num2istr(npks),extra=""
  804.     Mk(center,nrows);Mk(amplitude,nrows);Mk(width,nrows);Mk(area,nrows)
  805.     if(npks<1)
  806.         Abort "no peaks"
  807.     endif
  808.     DoWindow/F PeakFitGraph
  809.     if(V_Flag==0)
  810.         AppWv(g_w,g_wx,"");AppWv("W_PeakFit",wx,"");ColorWaves();DoWindow/C PeakFitGraph
  811.     else
  812.         NoCrs()
  813.     endif
  814.     
  815.     DoWindow/F PeakReportTable
  816.     if(V_Flag==0)
  817.         PeakReportTable()
  818.     endif
  819.  
  820.     | Separate the coefficients into separate waves
  821.     Variable pkTyp=$pdcw[5]
  822.     if(pkTyp==7) | MyPeak
  823.         sz=numpnts(W_MyPkCoefs);extra=""; |MyPeak
  824.         Mk(tw2,sz)
  825.         iterate (npks)
  826.             $tw2=$pcw[st+p]
  827.             n=MyCoefsToPeakSizes($tw2)
  828.             $center[i]=V_ctr;$amplitude[i]=V_amp;$width[i]=V_fwhm;$area[i]=V_area
  829.             st += sz
  830.         loop
  831.         text+=" "+S_PkT+" Peaks" | MyPeak
  832.     else
  833.         Abort "Unsupported Peak type"
  834.     endif
  835.     String cmd="Sort "
  836.     if(order==2)
  837.         cmd+="/R  "
  838.     endif
  839.     cmd+="$"+srt+",$center,$amplitude,$width,$area"+extra
  840.     Execute cmd        | Sort...
  841.     
  842.     DoWindow/F PeakReportLayout
  843.     if(V_Flag==0)
  844.         PeakReportLayout()
  845.     endif
  846.     if((basInf==1)%&(terms>0)%&(btyp>=0))
  847.         text+="\r\rBaseline function: "+S_funcs[btyp*10+1,btyp*10+7]
  848.         text+="\r\tK0 = "+num2str($pcw[1])
  849.         iterate (terms)
  850.             text+="\r\tK"+num2istr(i+1)+"= "+num2str($pcw[2+i])
  851.         loop
  852.     endif
  853.     ReplaceText/N=info text
  854.     ReplaceText/N=title "\JC"+title
  855.     Modify rows(PeakReportTable)=nrows
  856. End
  857.  
  858. Macro ShowFittedPeaks(w,wx,wb,pks,anno)
  859.     String w=g_w,wx=g_wx,wb=afp_b
  860.     Variable pks=afp_pks,anno=afp_anno
  861.     Prompt w,p_w,popup, WaveList("*",";","")
  862.     Prompt wx,p_wx,popup,calc+WaveList("*X", ";","")
  863.     Prompt wb,p_b,popup, none+"_From Fit Peaks_;"+WaveList("*base*",";","")
  864.     Prompt pks,"peak waves",popup,"off;on"
  865.     Prompt anno,"peak annotation tags",popup, "off;on"
  866.  
  867.     Silent 1;PauseUpdate
  868.     g_w=w;g_wx=wx;afp_anno=anno;afp_pks=pks;SBs(wb)
  869.     ChkLen(w,wx);ChkLen(w,wb)
  870.     String tdcw="W_PM",tcw="W_tmp",pkw,text,typ,tw2="W_tmp2"
  871.     if ((pks==1)%&(anno==1))
  872.         DoAlert 0,"Both tags and peaks are off!"
  873.     endif
  874.     if(numtype($pdcw[1])!=0)
  875.         Abort "Run FitPeaks first!"
  876.     endif
  877.     Variable s=$pdcw[1],en=$pdcw[2],terms,sz
  878.     sz=numpnts(W_MyPkCoefs)
  879.     Dup(pdcw,tdcw);Dup(pcw,tcw)
  880.     if(cmpstr(wb,"_From Fit Peaks_") == 0)
  881.         terms=$pdcw[4]
  882.         Redimension/N=(5+2) $tdcw;$tdcw[0]=2
  883.     else
  884.         $tdcw={2,s,en,0,0,NaN,sz}
  885.         $tcw={NaN,0}
  886.         terms = 0
  887.     endif
  888.     Variable dsz=2,st=2+terms,cst=2+$pdcw[4],dcst=5,n,y,npks=$pdcw[0]-1
  889.     Redimension/N=(st+sz) $tcw
  890.     iterate (npks)
  891.         $tcw[st,st+sz-1]=$pcw[cst+p-st]
  892.         $tdcw[5,6]=$pdcw[dcst+p-5]
  893.          if (pks==2)
  894.             NewWv(w,"_fpk");pkw=S_Wave;$pkw=NaN
  895.             if(exists(wx)==1)
  896.                 $pkw[s,en]=PolyMorph($tcw,$wx[p])
  897.             else
  898.                 $pkw[s,en]=PolyMorph($tcw,x)
  899.             endif
  900.             if(exists(wb)==1)
  901.                 $pkw[s,en]+=$wb[p]
  902.             endif
  903.             AppWv(pkw,wx,"");Modify mode($pkw)=1
  904.         endif
  905.         if($pdcw[5]==7) |MyPeak
  906.             typ="\JCVoigt\r"
  907.             Mk(tw2,sz)
  908.             $tw2=$tcw[st+p]
  909.             n=MyCoefsToPeakSizes($tw2)
  910.             y=$tcw[st+3]
  911.         else
  912.             Abort "Unsupported peak type"
  913.         endif
  914.         sprintf text, "\JL\Z09amp=%g\rcenter=%g\rwidth(+/-1%%)=%g\rarea=%g\ry=%g",V_amp,V_ctr,V_fwhm,V_area,y
  915.          if(anno==2)
  916.             FindPX(w,wx,V_ctr)
  917.             Tag/C/N=$("peak"+num2istr(g_ptgn))/A=MC  $w, V_peakX, typ+text
  918.             g_ptgn+=1
  919.         endif
  920.         dcst+=dsz;cst+=sz
  921.     loop
  922. End
  923.  
  924. Macro HideFittedPeaks()
  925.     Silent 1;PauseUpdate
  926.     String w,ws=WaveList("*_fpk",";","WIN:"+WinName(0,1))
  927.     Variable pos
  928.     do
  929.         pos=strsearch(ws,";",0)
  930.         if(pos<0)
  931.             break
  932.         endif
  933.         w=ws[0,pos-1]
  934.         ws[0,pos]=""
  935.         RmWv(w);KillWv(w)
  936.     while (1)
  937.     do
  938.         Tag/K/N=$("peak"+num2istr(g_ptgn))
  939.         g_ptgn-=1
  940.     while (g_ptgn>=0)
  941.     g_ptgn=0
  942.     ColorWaves()
  943. End
  944.  
  945. Proc FindPX(w,wx,xx)
  946.     String w,wx    | wx monotonic
  947.     Variable xx
  948.  
  949.     if(exists(wx)==1)
  950.         FindLevel/Q $wx, xx 
  951.         V_peakP = x2pnt($wx,V_levelX)
  952.         if(V_Flag==1)
  953.             if((xx<$wx[1])%&($wx[0]<$wx[1]))
  954.                 V_peakP=0
  955.             else
  956.                 V_peakP=numpnts($w)-1
  957.             endif        
  958.         endif
  959.         V_peakX = pnt2x($w,V_peakP)
  960.     else
  961.         V_peakP = x2pnt($w,xx)
  962.         V_peakX = xx
  963.     endif
  964. End
  965.  
  966. Proc ChkLen(w1,w2)
  967.     String w1,w2    | Wave names
  968.     if((exists(w1)==1) %& (exists(w2)==1))
  969.         if(numpnts($w1) != numpnts($w2))
  970.             Abort w1+" and "+w2+" have different number of points!"
  971.         endif
  972.     endif
  973. End
  974.  
  975. Proc Mk(w,n)
  976.     String w
  977.     Variable n
  978.     if(exists(w)==1)
  979.         Redimension/D/N=(n) $w;$w=NaN
  980.     else
  981.         Make/D/N=(n) $w=NaN
  982.     endif    
  983. End
  984.  
  985. Proc Dup(t,w)
  986.     String t,w
  987.     String n=""
  988.     if (exists(w)==1)
  989.         n=note($w) | keep w's note
  990.     endif
  991.     Duplicate/O $t,$w;Note/K $w
  992.     if( strlen(n)>0)
  993.         Note $w,n
  994.     endif
  995. End    
  996.  
  997. Proc SameLen(t,w)
  998.     String t,w
  999.     if(exists(w)!=1)
  1000.         Duplicate/O $t,$w
  1001.     else
  1002.         Redimension/D/N=(numpnts($t)) $w;CopyScales $t,$w
  1003.     endif
  1004. End    
  1005.  
  1006. Proc NewWv(srcw,sfx)
  1007.     String srcw,sfx
  1008.  
  1009.     Variable ii=1
  1010.     S_Wave=srcw[0,17-strlen(sfx)]+sfx
  1011.     if(exists(S_Wave)==1)
  1012.         String shw=S_Wave
  1013.         do
  1014.             S_Wave=srcw[0,14-strlen(sfx)]+num2istr(ii)+sfx
  1015.             ii+=1
  1016.         while ((exists(S_Wave)==1)*(ii<999))
  1017.     endif
  1018.     if(exists(srcw)==1)
  1019.         Duplicate/O $srcw,$S_Wave
  1020.     else
  1021.         Make/D/O $S_Wave
  1022.     endif
  1023.     Print "Created wave "+S_Wave
  1024. End
  1025.  
  1026. Proc RmWv(w)
  1027.     String w
  1028.  
  1029.     CheckDisplayed/W=$WinName(0,1) $w
  1030.     if(V_Flag==1)
  1031.         DoWindow/F $WinName(0,1);Remove $w
  1032.     endif
  1033. End
  1034.  
  1035. Proc KillWv(w)
  1036.     String w
  1037.     if(exists(w)==1)
  1038.         CheckDisplayed/A $w
  1039.         if(V_Flag==0)
  1040.             KillWaves $w
  1041.         endif
  1042.     endif
  1043. End
  1044.  
  1045. Proc AppCrs(w)
  1046.     String w
  1047.     ShowInfo;Cursor/P A,$w,0;Cursor/P B,$w,numpnts($w)-1
  1048. End
  1049. Proc NoCrs()
  1050.     HideInfo;Cursor/K A;Cursor/K B
  1051. End
  1052. Proc AppWv(w,wx,ax)
  1053.     String w,wx,ax
  1054.     
  1055.     String wn=WinName(0,1)
  1056.     if(strlen(wn)==0)
  1057.         ax = "Display $w"
  1058.     else
  1059.         DoWindow/F $wn
  1060.         ax="Append"+ax+" $w "
  1061.     endif
  1062.     if(exists(w)==1)
  1063.         CheckDisplayed $w
  1064.         if(V_Flag==0)
  1065.             if(exists(wx)==1)
  1066.                 ax+=" vs $wx"
  1067.             endif
  1068.             Execute ax
  1069.             ColorWaves()
  1070.         endif
  1071.     endif
  1072. End
  1073. Proc SBs(s)
  1074.     String s
  1075.     if((exists(s)==1)%|(cmpstr(s+";",none)==0))
  1076.         g_b=s;rb_b=s;fpks_b=s;afp_b=s
  1077.     else
  1078.     if(strsearch(S_funcs,s,0)>=0)
  1079.         fpks_b=s
  1080.     else
  1081.     if(cmpstr(s,"_From Fit Peaks_")==0)
  1082.         rb_b=s;afp_b=s
  1083.     else
  1084.     if(cmpstr(s,"_From Fit Baseline_")==0)
  1085.         fpks_b=s
  1086.     endif
  1087.     endif
  1088.     endif
  1089.     endif
  1090. End
  1091. Macro ShowFitResidual(w,wx,wfit,ow,extent,anno)
  1092.     String w=g_w,wx=g_wx,wfit=ar_wfit,ow=ar_ow
  1093.     Variable anno=ar_anno,extent=ar_ex
  1094.     Prompt w,p_wd,popup, WaveList("*",";","")
  1095.     Prompt wx,p_wxd,popup,calc+WaveList("*X", ";","")
  1096.     Prompt wfit,"fit wave",popup, WaveList("*Fit",";","")
  1097.     Prompt ow,"output residual wave",popup,new+WaveList("*",";","")
  1098.     Prompt extent,"extent of residual wave to measure",popup,"entire wave;between cursors;"+breg
  1099.     Prompt anno,"standard deviation lines annotations",popup,"± 1; ± 2;"+none
  1100.  
  1101.     Silent 1;PauseUpdate
  1102.     String wf=wfit[2,5]
  1103.     if(exists(ow)!=1)
  1104.         NewWv(w,wf+"Res");ow=S_Wave
  1105.     else
  1106.         SameLen(w,ow)
  1107.     endif
  1108.     g_w=w;g_wx=wx;ar_wfit=wfit;ar_ow=ow;ar_ex=extent;ar_anno=anno
  1109.     ChkLen(w,wx);ChkLen(w,wfit)
  1110.     Variable s=0,en=numpnts($w)-1
  1111.     if(extent==2) | use cursors
  1112.         CheckTwoCursors(w)
  1113.         s=V_start
  1114.         en =V_theEnd
  1115.     endif
  1116.     if(extent==3)
  1117.         if(exists(breg)!=1)
  1118.             Abort breg+" doesn't exist!"
  1119.         endif
  1120.         ChkLen(w,breg)
  1121.         $ow=($w[p]-$wfit[p]) / (numtype($breg[p])==0)
  1122.     else
  1123.         $ow=NaN;$ow[s,en]=$w[p]-$wfit[p]
  1124.     endif
  1125.     
  1126.     WaveStats/R=[s,en] $ow
  1127.     String rx="W_Resid"+wf+"X",ry="W_Resid"+wf+"Y"
  1128.     if(anno<3)
  1129.         Mk(rx,6);Mk(ry,6)
  1130.         if(exists(wx)==1)
  1131.             $rx[0,;4]=$wx[s]
  1132.             $rx[1,;4]=$wx[en]
  1133.         else
  1134.             $rx[0,;4]=pnt2x($w,s)
  1135.             $rx[1,;4]=pnt2x($w,en)
  1136.         endif
  1137.         $ry[0,;4]=(trunc(p/2)-1)*V_sdev*anno
  1138.         $ry[1,;4]=$ry[p-1]
  1139.         AppWv(ry,rx,"/R")
  1140.         Variable maxv=max(V_max,anno*V_sdev),minv=min(V_min,-anno*V_sdev)
  1141.         Modify zero(right)=1,minor(right)=1
  1142.         SetAxis right minv-3*(maxv-minv),maxv
  1143.     endif
  1144.     AppWv(ow,wx,"/R");Modify mode($ow)=3,marker($ow)=8,msize($ow)=1
  1145.     Label right "Residual (fit-data)"
  1146. End
  1147.  
  1148. Macro ShowOnlyDataAndBase()
  1149.     Silent 1;PauseUpdate
  1150.     String keep=g_w+";"+g_wx+";"
  1151.     keep+=g_b+";"+WaveList("*X",";","WIN:"+WinName(0,1))
  1152.     keep+=g_keep+";";g_keep=""
  1153.     String w,wn=WinName(0,1),ws=WaveList("*",";","WIN:"+wn)
  1154.     Variable pos,flg=0
  1155.     if(strlen(wn)!=0)
  1156.         do
  1157.             pos=strsearch(ws,";",0)
  1158.             if(pos<0)
  1159.                 break
  1160.             endif
  1161.             w=ws[0,pos-1]
  1162.             ws[0,pos]=""
  1163.             if(strsearch(keep,w+";",0)<0)
  1164.                 RmWv(w);flg=1
  1165.             endif
  1166.         while (1)
  1167.         if(flg)
  1168.             ColorWaves()
  1169.         endif
  1170.     endif
  1171. End
  1172.  
  1173. Macro ColorWaves() : GraphStyle
  1174.     Silent 1;PauseUpdate
  1175.     Modify/Z lSmooth=1,lHair=0.5,minor(bottom)=1
  1176.     Modify/Z rgb[0]=(65535,0,0)
  1177.     Modify/Z rgb[1]=(0,0,65535),rgb[2]=(3009,65535,1882),rgb[3]=(0,0,0)
  1178.     Modify/Z rgb[4]=(0,0,65535),rgb[5]=(63953,3661,65535)
  1179.     Modify/Z rgb[6]=(37510,1,1),rgb[7]=(27232,40528,22540),rgb[8]=(1531,1314,28456)
  1180.     Legend/C/N=default ""
  1181. End
  1182.  
  1183.  
  1184. Proc CheckTwoCursors(w)
  1185.     String w    | the y wave
  1186.     Variable/G V_start,V_theEnd,V_startX,V_theEndX | point indices, x vals
  1187.     String wa=CsrWave(A),wb=CsrWave(B),t="cursors on target graph "
  1188.     if((strlen(wa)==0)%|(strlen(wb)==0))
  1189.         Abort "expecting two "+t
  1190.     endif
  1191.     V_start=pcsr(A);V_theEnd=pcsr(B)
  1192.     V_startX=hcsr(A);V_theEndX=hcsr(B)
  1193.     Variable x0=pnt2x($w,0),x1=pnt2x($w,1),ok
  1194.     ok = (numpnts($wa)==numpnts($w)) %& (numpnts($wb)==numpnts($w)) 
  1195.     ok = ok %& (x0==pnt2x($wa,0)) %& (x1==pnt2x($wa,1))
  1196.     ok = ok %& (x0==pnt2x($wb,0)) %& (x1==pnt2x($wb,1)) 
  1197.     if(ok)
  1198.         if(V_start>V_theEnd)
  1199.             x0=V_Start
  1200.             V_start=V_theEnd
  1201.             V_theEnd=x0
  1202.             x0=V_startX
  1203.             V_startX=V_theEndX
  1204.             V_theEndX=x0
  1205.         endif
  1206.         if(V_start==V_theEnd)
  1207.             Abort t+"are too close together!"
  1208.         endif
  1209.     else
  1210.         Abort t+"must be on wave "+w+" or one like it"
  1211.     endif
  1212. End
  1213.  
  1214. Window DefinePeaks() : Table
  1215.     PauseUpdate; Silent 1        | building window...
  1216.     Edit W_MakeCentersX.y,W_MakeAmpsY.y,W_MakeWidthsX.y as "DefinePeaks"
  1217. EndMacro
  1218.  
  1219. Window PeakEstimates() : Table
  1220.     PauseUpdate; Silent 1        | building window...
  1221.     Edit W_EstCentersP.y,W_EstCentersX.y,W_EstAmpsY.y as "Peak Estimates"
  1222.     Append W_EstWidthsX.y,W_EstEdgesP.y
  1223.     Modify width(Point)=64
  1224. EndMacro
  1225.  
  1226. Window PeakReportTable() : Table
  1227.     PauseUpdate; Silent 1        | building window...
  1228.     Edit W_PkCenters.y,W_PkAmps.y,W_PkFWHMs.y,W_PkAreas.y
  1229. EndMacro
  1230.  
  1231. Window PeakReportLayout() : Layout
  1232.     PauseUpdate; Silent 1        | building window...
  1233.     Layout /W=(8,41,472,332) PeakReportTable(124,388,518,438)/O=2 as "Peak Report Layout"
  1234.     Textbox /N=title/A=LT/X=39.8182/Y=5.35714 "\JCYour Report Title Here"
  1235.     Modify left(title)=250,top(title)=71,width(title)=107,height(title)=6,frame(title)=4
  1236.     Textbox /N=info/A=LT/X=39.0909/Y=43.8187 "1 Hanning Peaks"
  1237.     Modify left(info)=246,top(info)=351,width(info)=75,height(info)=10,frame(info)=1
  1238.     Append PeakFitGraph(63,116,555,325)/O=1
  1239.     Modify mag=1
  1240. EndMacro
  1241.  
  1242.